Bending Modulus of Lipid Membranes from Density Correlation Functions

The bending modulus κ quantifies the elasticity of biological membranes in terms of the free energy cost of increasing the membrane corrugation. Molecular dynamics (MD) simulations provide a powerful approach to quantify κ by analyzing the thermal fluctuations of the lipid bilayer. However, existing methods require the identification and filtering of non-mesoscopic fluctuation modes. State of the art methods rely on identifying a smooth surface to describe the membrane shape. These methods introduce uncertainties in calculating κ since they rely on different criteria to select the relevant fluctuation modes. Here, we present a method to compute κ using molecular simulations. Our approach circumvents the need to define a mesoscopic surface or an orientation field for the lipid tails explicitly. The bending and tilt moduli can be extracted from the analysis of the density correlation function (DCF). The method introduced here builds on the Bedeaux and Weeks (BW) theory for the DCF of fluctuating interfaces and on the coupled undulatory (CU) mode introduced by us in previous work. We test the BW-DCF method by computing the elastic properties of lipid membranes with different system sizes (from 500 to 6000 lipid molecules) and using coarse-grained (for POPC and DPPC lipids) and fully atomistic models (for DPPC). Further, we quantify the impact of cholesterol on the bending modulus of DPPC bilayers. We compare our results with bending moduli obtained with X-ray diffraction data and different computer simulation methods.

The simulation contained 2304 DPPC and 2304 cholesterol molecules and were performed with the GROMACS simulation package 6 .

B. All-Atom DPPC
The All-atom simulations of the DPPC bilayers were performed using the CHARMM36 forcefield 7 at 323.15 K. We used the v-rescale thermostat with a time constant of 1 ps, and the Parrinello-Rahman barostat was employed to perform the simulations at constant surface tension equal to zero, with a pressure normal to the membrane plane of 1 bar. The coupling constant for the barostat was set to 5 ps. We employed the Particle Mesh Ewald method with a cutoff in real space of 1.2 nm, and the Lennard-Jones interactions were truncated at 1.2 nm using a switching function between 1 and 1.2 nm. The initial configurations were generated using CHARMM-GUI 8- 10 We performed simulations of bilayer containing 1000 and 1800 DPPC molecules.
In Table I we list the main parameters of the simulations.

II. ISM RESULTS FOR THE DPPC-CHOLESTEROL MEMBRANE
As the DPPC-Cholesterol membrane had not been previously analyzed with the ISM procedure, we include here the ISM results for that case. Figure 1 shows the mean density profile, the intrinsic profile and the mean density profile obtained via the convolution: between an intrinsic density profile, ρ I (z), and the probability distribution, P(ξ), of the local values of the intrinsic surface z = ξ( x). Dashed green line: mean density profile obtained from the Capillary Wave Theory (eq. 1 ), with the probability distribution P (ξ) corresponding to the MD-ISM used to get ρ in (z; q u ).
In Figure 2 we show the surface tension dependence with q for the tensionless DPPC-Cholesterol membrane. In addition to the modes shown in the Figure 1 of the main text we also include the peristaltic (P) mode 2 , which is given by the local half-thickness The peristaltic mode does not show the divergence |ξ x q | 2 ∼ 1/q 2 characteristic of fluctuating surfaces (x = m or U ), therefore it provides a more natural way of defining u P (q) = kT /( |ξ P q | 2 A o ), with a finite limit u P (0). To obtain a q-dependent surface tension for the P mode, similar to eq (5) of the main text, we define γ P (q) = u P (q)/q 2 and which is shown in Figure 2. red the peristaltic (uncoupled) mode γ P ef f (q); blue the undulatory mode γ U (q). and black the monolayer mode γ m (q). The dashed are a guide for the eyes. Figure 2 , the addition of cholesterol at (1:1) concentration in DPPC, expands the range of q-vectors ( γ CU (q) < γ P ef f (q) ) where the fluctuations of the two lipid layers are correlated.

As shown in
Of the three modes; U, CU and P, only two are independent and are related to each other by equation 2 : As explained in the main text, the function D(q) defined by the ratio between the U and CU surface tensions, is expected to behave as Within of the ISM formalism we can evaluate D(q) directly, and check that indeed it has no q 2 term in its Taylor expansion . In the next figure we show the ISM D(q) for the free DPPC membranes.   red the peristaltic (uncoupled) mode γ P ef f (q); blue the undulatory mode γ U (q). and black the monolayer mode γ m (q). The lines are a guide for the eyes.
Compared with the MARTINI forcefield model, the all-atom forcefield for DPPC gives a much narrower window ( γ CU (q) < γ P ef f (q) ) of q values where both lipid layers fluctuate together.

FUNCTIONS OF LARGE MEMBRANES
As shown in Figure 6, for large lateral sizes the mean density profile of the tensionless POPC membrane features two wide peaks that spread over a larger z region, overlapping with each other. The impact of the larger fluctuations in these large membranes can also be observed in the density correlation functions, with the four quadrants spreading over larger (z 1 ,z 2 ) regions and overlapping with each other (see Figure 7) . We have used the same scale for the axes (x,y) as in Figure 2 in the main article to facilitate comparison. We show in Figure 8 the impact of the surface tension. Upon applying a tension on the membrane the density correlation shows again four separate quadrants. Notice that the second column of Figure 8 for the membrane with 4000 POPC molecules per layer corresponds to the same q = 0.118nm −1 as the first column of Figure 2 of the main article, for a membrane with 1000 POPC molecules under the same tension. Wertheim's term (n = 1 in the BW series) gives the same (size independent) prediction for both systems, but the MD results for G(z, z , q) are clearly different and the full BW series describes well the influence membrane size on the DCF for a given q.

V. THE STRUCTURE FACTOR
In Figure 9 we compare the simulated structure factor S(q, q z ) with the BW structure factor S BW (q, q z ) for a POPC membrane under tension, with 1000 lipids per layer, at the lowest q = 2π/L x compatible with the simulation box. As we can see at this low wavevector, q = 2π/L x = 0.237nm −1 , the BW S BW (q, q z ) represents more accurately the simulated S(q, q z ) results than at q = 10π/L x = 1.18nm −1 shown in Figure 3 of the main text.
In Figure 10 we compare the interlayer simulation structure factor S +− (q, q z ) with the interlayer BW structure factor S +− BW (q, q z ) to q = 10π/L x for the DPPC-Cholesterol (top panel) and for the tensionless DPPC membranes (bottom panel). This figure is equivalent to Figure 4 of the main text, but now for the tensionless DPPC and DPPC-Cholesterol membranes. Therefore the BW predictions using γ U (q) and γ CU (q) are much more similar to each other than for a pure DPPC membrane (at the same q values). The Figure 12 is the same Figure 6 of the main text, i.e the inverse of the coupled undulatory surface tension, γ CU (q), for the DPPC membranes analyzed in this work, but now including the inverse of the usual undulatory surface tension, γ U (q), obtained with the ISM. The γ U (q) function is the usual way to obtain the bending but , as we discuss in the main text, it is not easily obtained from the density correlation function since it contains a large contribution the non-BW background. The BW-DCF route introduced in this work is only applicable to the coupled undulatory mode, γ CU (q).

VII. ROBUSTNESS OF THE FITS TO q 2 γ CU (q) FOR THE MARTINI FREE POPC MEMBRANES
In this section we check, for the MARTINI free POPC membranes, the robustness with the lateral size of the membranes. of the fits of q 2 γ CU (q) to the parametric form given by eq. (19) of the main text. For this we have used the simulations MARTINI free POPC membranes carried out in our previous work 3 . As shown in table II the agreement between the BW-DCF and ISM data is good across all the systems studied. We note that use of the BW-DCF procedure requires sufficiently large system sizes, typically N > 1500. For smaller system sizes the number the points within of the relevant fitting interval 0 < q < q d (where the CWT and the low-q approximation to D(q) (see eq. (4) do not apply) is reduced drastically. However, if we use the ISM values for D(q) instead of the low q approximation we can obtain acceptable estimation for the bending modulus using systems containing N = 500 lipids. This analysis falls outside the main objective of this work.